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INTRODUCTION 


Spall  is  a  widely  observed  phenomenon  accompanying  underground  nuclear 
explosions,  and  its  seismic  consequences  are  therefore  of  considerable  interest  in  the 
context  of  verification  of  test  limitation  treaties.  The  proper  representation  of  spall  is  of 
importance  for  forward  and  inverse  seismic  modeling  of  explosions  for  yield  estimation 
and  discrimination  studies. 

While  spall  is  expected  to  have  negligible  effect  on  long  period  seismic  waves  (see 
Day  et  al.,  1983),  several  studies  suggest  that  it  is  a  significant  contributor  to  the  short 
period  seismic  signal  from  explosions.  If  this  is  the  case,  spall  may  complicate  explosion 
yield  estimates,  or  at  least  influence  our  interpretation  of  empirically  based  yield  estimation 
formulas. 

For  example,  the  isotropic  seismic  moment  provides  one  estimate  of  explosion 
yield.  Patton  (1988)  used  regional,  higher  mode  Rayleigh  waves  in  the  0.2  to  0.5  Hz  band 
to  estimate  the  isotropic  moment  of  the  underground  explosion  HARZER.  He  concluded 
that  the  isotropic  moment  estimate  v/as  sensitive  to  assumptions  about  the  amount  of  spall 
and  its  efficiency  in  generating  seismic  waves  in  this  frequency  band. 

As  a  second  example,  several  studies  (e.g.,  McLaughlin,  et  al.,  1988;  Taylor  and 
Randall,  1989)  have  suggested  that  spall  could  be  a  significant  contributor  to  the  regional 
phase  Lg,  which  currently  shows  great  promise  as  a  low-variance  yield  estimator  for 
explosions  at  the  Soviet  Shagan  River  test  site  (Hansen  et  al.,  1990).  This  suggestion 
arises  because  spall  has  a  deviatoric  source  component,  which  may  significantly  enhance 
the  production  of  short  period  SV  waves,  compared  with  the  predictions  of  purely  isotropic 
explosion  source  models. 

Teleseismic  P  wave  amplitudes  provide  a  third  basis  for  yield  estimation.  At 
teleseismic  distances,  the  effect  of  spall  is  to  partially  cancel  the  pP  phase,  replacing  it  with 
an  attenuated  and  delayed  phase  which  is  not  a  replica  of  the  direct  P  wave  (e.g.,  Day  et  al., 
1986).  Schlittenhardt  (1990)  presents  theoretical  calculations  showing  that  spall  can 
potentially  influence  teleseismic  P  wave  amplitudes,  given  plausible,  but  uncertain, 
assumptions  about  the  characteristic  mass  and  momentum  of  spall. 
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Day  et  al.  (1983)  (hereafter  referred  to  as  paper  D83)  presented  a  theoretical 
argument  that  surface  waves  of  period  20  seconds  and  longer  are  not  significantly  affected 
by  spall.  These  authors  then  proposed  that  the  seismic  consequences  of  spall  be  modeled 
by  a  shallow,  horizontally  oriented  tension  crack.  They  analyzed  the  implications  of  this 
model  for  surface  waves  with  wavelength  long  compared  to  the  crack  dimension  and  depth, 
and  showed  that,  in  this  approximation,  the  tension  crack  model  was  equivalent  to  a 
vertically  oriented  point  force  acting  at  the  free  surface,  with  the  force  time  history 
proportional  to  the  crack-separation  acceleration  history. 

The  resulting  point  force  model  has  the  property  that  the  model  parameters  are,  at 
least  in  principle,  observable  quantities.  The  spall  mass  and  momentum,  for  example,  can 
be  estimated  from  near-field  strong  motion  recordings  at  the  U.S.  Nuclear  Test  Site  (NTS) 
(e.g.,  Viecelli,  1973;  Sobel,  1978;  Patton,  1990).  This  property  of  relating  seismic  waves 
to  near-field  observables  has  proven  attractive,  and  the  point  force  model  of  D83  has  been 
applied  extensively  (e.  g.,  Stump,  1985;  Patton,  1988;  Taylor  and  Randall,  1989), 
sometimes  under  conditions  which  depart  significantly  from  those  for  which  it  was 
originally  conceived,  namely,  low-frequency  surface  waves. 

For  this  reason,  we  conclude  that  a  useful  purpose  is  served  by  presenting  a 
derivation  of  the  point  force  model  which  is  more  general  than  the  original  derivation.  The 
new  derivation  brings  out  more  clearly  the  relation  of  the  point  force  representation  to  the 
tension-crack  model,  from  which  it  arises  by  approximation.  The  derivation  also  clarifies 
the  relationship  between  point  force  and  moment  tensor  representations  for  spall,  a  result 
which  has  implications  for  the  explosion-source  inverse  problem.  We  then  examine 
numerically  the  limitations  of  the  point-force  representation.  Finally,  we  show  that  some  of 
the  assumptions  inherent  in  the  original  tension-crack  model  can  be  relaxed  without 
requiring  revision  of  the  point-force  and  moment  tensor  representations. 

This  paper  addresses  the  seismic  representation  of  spall  primarily  within  the 
confines  of  the  point  source  approximation,  i.e.,  wavelength  long  compared  with  the 
characteristic  source  dimensions.  This  is  the  framework  within  which  most  previous 
forward  modeling  and  source  inversion  has  been  conducted.  The  important  ancillary 
effects  of  lateral  finiteness  of  the  spall  source  will  be  discussed  in  a  subsequent  paper 
(Barker  and  Day,  in  preparation). 
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TENSION  CRACK  MODEL 


D83  proposed  that  spall  be  represented  as  a  horizontally  oriented  tension  crack  that 
opens  and  closes  in  the  vertical  direction  (see  Figure  1).  With  the  aid  of  a  Green's  tensor 
G,  the  displacement  field  u  due  to  a  tension  crack  can  be  written  as  a  surface  integral  over 
the  crack  surface  £  (AM  and  Richards,  1980,  p.  39): 

Mi  {x)=  J  VjiO  SUk{£)  Cjkpqi £)  G/M(*,$  d£  , 

(1) 

where  v  is  the  unit  normal  to  the  crack,  £  is  the  general  position  on  £,  8u  is  the  spall 
separation, 

i.  e.,  the  displacement  discontinuity  across  the  crack,  and  C  is  the  elastic  tensor.  For 
horizontal  crack  orientation,  vertical  spall  separation,  and  crack  depth  h,  (1)  can  be  written 
in  the  form 


m«  (*)  = 


l 


Mpq(V)Gip,q(X‘,Tf)  d^T]  , 


(2) 


where  i]  is  the  general  position  in  the  source  volume  V,  and  m  is  the  volumetric  moment 


tensor  density, 


mpq(V)  =  8u3(Vuri2)C33pq  %V3~h) 


( S(  ■)  is  the  Dirac  delta  function).  Then,  in  the  point  source  approximation  (wavelength 
large  relative  to  the  maximum  dimension  of  £),  (2)  reduces  to 


Mi  (x)  =  MpqGip>q(x;0,0,h)t  (4) 

r 

where  the  moment  tensor  M  (volume  integral  of  m)  for  an  isotropic  earth  model  has  the  ~ 
following  matrix  components: 
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M  =  &73  A  0  A.  0 
0  0  A  +  2/i 

Here  5uz  is  the  spall  separation  averaged  over  the  crack  area  A,  and  X  and  (i  are  the  Lame 
constants.  It  will  be  convenient  to  refer  to  (4)  and  (5)  as  the  spall  moment  tensor 
representation. 


POINT  FORCE  REPRESENTATION 
An  expression  equivalent  to  (2)  is 

Mi  (X)  =  Su3(T}hT]2)  r33(7h,772,fc)  dT]Xdt]2  , 

JJ  (6) 

in  which  we  introduce  the  notation 

733(7?)  =  CyipqG  ip^q{x\7f)  (7) 

This  notation  is  introduced  to  emphasize  that  this  factor  can  be  interpreted,  by  reciprocity, 
as  a  component  of  the  stress  tensor  T  induced  at  the  (source)  location  7)  in  response  to  a 
point  force  (in  the  i  direction)  acting  at  the  (receiver)  location  x.  That  is,  T  is  simply  the 
stress  tensor  derived  from  the  elastodynamic  displacement  field  of  a  point  force  (note  that 
this  usage  of  the  symbol  T  differs  from  that  of  chapters  2  and  3  of  Aki  and  Richards,  in 
which  T  represents  a  general  traction  vector).  We  expand  the  ?)3  dependence  of  T33  in  a 
Taylor  series  about  the  free  surface,  7]3  =  0: 

733(f?i, t?2»h)  =  733(771,772,0)  +  hr33,3(  771,772,0)  +  0(h2)  (3) 

The  first  term  is  zero  by  virtue  of  the  free  surface  boundary  condition.  Furthermore,  T 
satisfies  the  equation  of  motion,  so  (in  any  region  which  excludes  the  point  x)  we  have 

TsjM  =  -0)2pGn(x',Ti)'(9) 

The  free  surface  boundary  conditions  ensure  that  TV  vanishes,  for  all  j  and  for  k=l,2,  in 
the  limit  that  7j3  goes  to  zero.  Hence,  the  limiting  form  of  the  momentum  equation  (9)  is 
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733,3(^1^2,0)  =  -©2pG!i3(x;7J1,772,0)  , 


(10) 


« 


& 


Substituting  (10)  into  the  Taylor  series  (8),  and  the  result  into  the  surface  integral 
representation  (6),  yields,  to  first  order  in  the  spall  depth  h, 


Ui  (x)  =  -h 


pcfi8u3(7h,ri2)  Ga(x;rh,ii2, 0)  dTjidi]2 


(11) 


In  (11),  the  Green's  function  is  evaluated  for  source  points  at  the  free  surface,  so 
that  its  coefficient 


gs  =  -h  p  ofi8u2(rih7]2)  (12) 

can  be  interpreted  as  a  normal  traction  applied  at  the  earth's  surface  (not  the  spall  surface), 
with  amplitude  given  by  the  areal  density  of  the  spall  layer,  ph,  and  time  history 
proportional  to  the  acceleration  of  the  spall  layer  relative  to  the  substrate.  The  integral  (11) 
thus  represents  spall  as  a  distributed  vertical  surface  traction;  this  is  an  accurate  alternative 
to  the  moment  tensor  density  representation  (2)  as  long  as  the  higher  order  terms  in  the 
Taylor  series  (8)  are  negligible,  i.e.,  when  the  wavelength  is  long  compared  to  the  source 
depth.  This  criterion  will  frequently  be  satisfied  in  practice;  for  example,  for  a  depth  of  150 
meters  and  wavespeed  of  3000  m/sec,  a  quarter  wavelength  criterion  for  spall  depth 
corresponds  to  5  Hz  seismic  waves.  Numerical  experiments  discussed  below  indicate  that, 
in  practice,  (1 1)  is  actually  accurate  for  spall  depths  up  to  about  half  a  wavelength. 

Finally,  in  the  point  source  approximation,  (1 1)  becomes  what  we  will  call  the  spall 
point  force  representation: 


m  (x)  =  Fs  Ga(x; 0),  (13) 

where  the  point  force  amplitude  Fs  is  the  surface  integral  of  as,  and  equals  the  product  of 
the  spall  mass  ms  and  the  crack  separation  acceleration : 


Fs  =  -ms  ofi  8u3  t 

(14) 

tns=  phA 

(15) 
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Strong  motion  observations  show  that  the  lateral  extent  of  spall  typically  exceeds  its 
depth  by  a  substantial  amount.  For  example,  Patton  (1990)  finds  that,  on  average,  spall 
radius  is  roughly  a  factor  of  6  greater  than  spall  depth  at  the  Pahute  Mesa  area  of  NTS.  As 
a  consequence,  lateral  finiteness  effects  will  vitiate  the  point  source  approximation  long 
before  the  half  wavelength  criterion  for  the  spall  depth  is  violated.  Thus,  the  spall  moment 
tensor  and  spall  point  force  are,  in  practice,  equivalent  mathematical  representations  of  the 
tension  crack  model.  From  (5),  (14),  and  (15),  it  follows  that  the  source  spectra  M  and  Fs 
are  related  by 


M(a»  =  --a-FV(s>) 
hor 


*/(X+2  n)  0  0 

0  V(A  +  2 pL)  0 

0  0  1 


(16) 


where  a  is  the  P  wave  speed.  The  spall  moment  tensor  time  history  is  cfi/h  times  the 
doubly  time  integrated  point  force  time  history. 

NUMERICAL  EXAMPLE 


We  use  regional-distance  synthetic  seismograms  to  verify  the  assertion  that  the 
point  force  and  moment  tensor  formulations  are.  equivalent  throughout  the  frequency  range 
in  which  lateral  finiteness  can.  be  neglected.  Our  focus  on  regional  distance  reflects  the 
importance  of  that  distance  range  for  the  seismic  verification  problem.  Moreover,  regional 
seismograms  contain  a  complex  mixture  of  wave  types,  and  provide  an  appropriate  test  of 
generality  of  the  purported  equivalence  of  point  force  and  moment  tensor.  Synthetics  are 
computed  using  the  PROSE  code  developed  by  Apsel  (1979). 

Figure  2  compares  synthetic  vertical  component  displacements  at  300  km  range  for 
the  two  source  representations.  The  spall  separation  time  history  is  a  delta  function,  the 
earth  model  (Table  1)  is  the  Eastern  Kazakhstan  crustal  model  of  McLaughlin  et  al.  (1988), 
and  the  synthetics  are  complete  for  the  0-5  Hz  frequency  band.  The  point  force  and 
moment  tensor  representations  give  indistinguishable  results  for  all  seismic  phases  over  this 
pass  band  when  the  spall  depth  is  100  meters.  For  larger  spall  depths,  the  agreement  at 
high  frequency  only  gradually  degrades. 
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Figure  3  compares  the  Fourier  spectra  for  the  two  representations,  for  windows 
taken  around  the  Pg,  Lg,  and  Rg  phases,  respectively.  Evidently,  in  the  absence  of  lateral 
finiteness  effects,  the  spall  point  force  representation  is  an  excellent  approximation  to  the 
spall  moment  tensor  representation  for  spall  depths  up  to  at  least  one  half  wavelength. 
Patton  (1990)  found,  for  Pahute  Mesa  explosions  below  the  water  table,  that  maximum 
spall  depth  was  approximately  60  m/kt1^3.  For  a  125  kt  event,  for  example,  the  point  force 
representation  would  not  be  expected  to  deviate  significantly  from  the  moment  tensor 
representation  for  frequencies  less  than  roughly  2.5  Hz  (assuming  a  near-surface  S  wave 
velocity  of  about  1.6  km/sec,  following  Bache  et  al.,  1978).  As  already  noted,  the  point 
source  approximation  is  invalidated  by  lateral  fmiteness  effects  at  considerably  lower 
frequencies  than  this. 


GENERALIZED  TENSION  CRACKS 

In  the  tension  crack  model  which  is  the  basis  of  the  foregoing  analysis,  material 
failure  is  confined  to  a  horizontal  plane.  It  is  possible,  however,  that  in  actual  explosions, 
the  depth  of  the  detachment  horizon  may  vary  laterally  across  the  spall  region.  For 
example,  Stump  (1985)  found  evidence  of  variable-depth  spall  in  sub-surface  accelerometer 
recordings  of  a  buried  chemical  explosion.  Furthermore,  some  numerical  simulations  of 
buried  explosions  indicate  spalling  over  a  roughly  conical  surface,  deepest  beneath  ground 
zero  and  shallowing  with  increasing  radial  distance  from  ground  zero  (see,  e.g.,  Figure  1 
of  Walton  and  Heuze,  1989).  In  such  cases,  vertical  movement  of  the  spall  mass  implies 
both  nonnal  and  tangential  components  of  relative  motion  across  the  spall  surface.  It  might 
be  supposed  that  such  an  effect  would  require  modification  of  the  spall  representations 
discussed  in  the  previous  sections. 

However,  a  simple  analysis  shows  that  the  spall  moment  tensor  and  spall  point 
force  representations  derived  from  the  horizontal  tension  crack  model  will  usually  be  valid 
even  when  the  spall  surface  2  is  non-planar,  provided  (i)  spall  separation  is  predominantly 
vertical,  and  (ii)  the  mean  dip  of  2  is  small  (but  arbitrarily  rough  relief  is  permitted  on  2). 
Assumption  one  is  supported  by  strong  motion  observations  in  the  spall  zone.  Likewise, 
assumption  two  is  supported  by  the  observation  that  lateral  extent  of  spall  generally  exceeds 
its  depth  by  a  substantial  factor  (recall  that  Patton,  1990,  infers  a  radius  to  depth  ratio  of 
roughly  6).  Furthermore,  any  axisymmetric  spall  surface  has  zero  mean  dip,  and  any  spall 
surface  which  terminates  by  intersecting  the  earth’s  surface  has,  neglecting  asymmetric 
topography,  zero  mean  dip. 
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To,  verify  the  assertion  that  the  moment  tensor  and  point  force  derived  from  the 
horizontal  tension  crack  have  this  more  general  validity,  we  take  the  spall  depth  to  be  a 
given  function  h(T]i,Tj2)  of  the  horizontal  coordinates,  so  that  the  spall  surface  is  specified 

by 


Tfe  -  h(rii,ri2)  =  0  f 

(17) 

The  representation  theorem  (1)  becomes 

Wj  (x)=  J  8ll3  Yf  Cjipq  G'ip,q 

(( u  -  r  dlhdll2  • 

-  I  I  d«3  Vj  ^j2pq  GiPtq  ^ 

(18) 

(the  double  integration  is  over  the  projection  of  £  onto  the  771,  %  plane)-  We  decompose 
the  summation  on  j ,  which  is  implicit  in  (18),  into  two  parts.  The  first  part,  Iv, 
corresponds  to  j  =  3  (vertical  component  of  the  normal  vector  v),  the  second  part,  ///( 
corresponds  to;  =  1, 2  (horizontal  components  of  v): 

Ui(x)=Iy  +  In , 

(19) 

in  which 

Iv  =  JJ  Su3  733(771,%, KVbmMrhdfy 

(20) 

U{  -  Jj  81I3  T3a~dT]ldT]2 

(21) 

where  a  =  1,2,  and  the  summation  convention  on  a  applies.  The  first  term,  Iy ,  is 
analogous  to  our  result  for  the  simple  tension  crack,  and  can  be  transformed  to  a  surface 
traction  representation  by  the  same  steps  leading  to  (1 1-12): 


/v=  oi(7h,7?2)  Ga(x;-rjh 772,0)  drjidm 

JJ  ,  (22) 
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where  (Js  is  the  surface  density  of  the  (variable-depth)  spall  layer  times  the  spall 
acceleration, 


=  -pKVi,V2)co28u3(ih,i]2) 


In  the  point  source  approximation,  we  obtain  the  analogue  of  (13-15), 


Iv  =  Fs  Gi 3(*;0) } 


where 


=  ll  <?s(rib 


m)  dJhdm 


The  second  tenn,  In  (Equation  21),  can  be  transformed  using  the  identity 

va_  dh 
V3  dVa 


to  give 


■// 


///  =  |  I  Su2  T3odah  dmdm 


Since  T^a  vanishes  at  h  =  0,  (27)  is,  to  first  order  in  h, 


=  jj  da</>( 


h  8u3  T3a, 3  drj\drj2 


<P  8113  73a, 3  driidri2 


where 


0  =  1*2 
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In  the  point  source  approximation,  (28)  becomes 


■II 


Ih  =  730,3  do0  S113  dl}\dTJ2 


(30) 


and  if  we  add  the  assumption  that  SU3  is  uncorrelated  with  da  <p,  then  we  can  factor  the 
mean  spall  separation  SI3  out  of  the  integral  in  (30).  The  result  is 


/;/=  ^<5«3T3 


h 2  Ya  dl 

> 


(31) 


where  the  line  integral  is  over  the  periphery  of  the  projection  of  the  spall  surface  £  onto  the 

A 

horizontal  plane,  and  7  is  the  unit  (two-dimensional)  outward  normal  to  this  curve  (lying  in 
the  7/i,  772  plane)-  The  line  integral  in  (31)  vanishes  whenever  the  spall  surface  is  bounded 
by  a  level  curve,  and  will  be  small  whenever  the  mean  dip  of  the  spall  surface  is  small.  The 
latter  interpretation  follows  if  h  is  expanded  to  first  order  about  the  mean  depth,  ho.  Then 
the  integral  reduces  to  hoAda  ,  where  A  is  the  (projected)  area  of  spall,  and  da  is  the 
average  value  of  the  slope  dah.  Thus,  it  is  likely  that  In  can  be  neglected  in  most  cases  of 
interest.  The  remaining  term,  /y,  is  identical  to  the  point  force  representation  for  the  simple 
crack,  as  shown  by  comparing  (24-25)  with  (13-15)  (although  we  can  no  longer  factor  Fs , 
Equation  25,  into  the  product  of  spall  mass  and  average  acceleration,  unless  and  h  are 
uncorrelated  in  771, 772). 


DISCUSSION 

D83  derived  the  point  force  representation  (Equations  13  - 15)  for  the  special  case 
of  Rayleigh  waves.  The  present  derivation  justifies  the  subsequent  application  of  the  point 
force  representation,  by  numerous  authors,  to  model  other  seismic  phases  (i.e.,  the 
approach  is  justified  to  the  extent  that  the  generalized  tension  crack  approximates  the 
geometry  uf  actual  spall,  the  separation  function  Suj  is  appropriately  prescribed  to  represent 
the  kinematics  of  spall,  and  the  wavelength  limitation  is  honored).  The  derivation  also 
serves  to  underscore  that  fact  that,  while  the  spall  moment  tensor  is  located  at  spall  depth, 
the  spall  point  force  is  located  at  the  earth's  surface. 
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The  derivation  makes  clear  the  redundancy  of  representing  spall  using  both  the  spall 
moment  tensor  and  the  spall  point  force  concurrently;  since  the  two  are  equivalent,  any 
attempt  to  infer  simultaneously  the  spectra  of  both  of  these  source  terms  from  seismic 
waveforms  will  lead  to  a  singular  inverse  problem.  This  was  previously  noted  by  Stump 
(1987,  1990)  on  the  basis  of  numerical  experiments.  Generally,  a  singular  inverse 
problem  will  result  if  the  point  force  representation  is  combined  with  any  moment  tensor 
representation  from  which  (5)  can  be  constructed  by  linear  combination. 

Equation  16,  which  expresses  the  correspondence  between  the  source  spectra  M(ai) 
and  Fs{cd),  has  an  additional  consequence  which  should  be  noted,  If  M(co)  is  estimated  by 
inversion  of  seismic  data  (e.  g.,  Stump,  1988,  1990;  Johnson,  1988),  the  numerical  values 
associated  with  M(co)  are  expected  to  depend  inversely  upon  the  source  depth  which  is 
assumed  for  the  inversion  (whereas  inversion  for  Fs(co)  does  not  incorporate  an  assumption 
of  source  depth).  This  point  becomes  important  when  inversions  based  on  the  two 
alternative  representations  are  compared  quantitatively ,  as  they  sometimes  have  been  in  the 
seismic  literature  on  spall. 
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FIGURE  CAPTIONS 


Figure  1.  Geometry  of  the  tension  crack  spall  model. 

Figure  2.  Synthetic  seismograms  for  the  point  force  and  moment  tensor  representations. 
Sources  for  the  4  depths  h  have  been  scaled  to  constant  spall  mass  phA  . 

Figure  3.  Fourier  spectral  ratios  for  the  synthetics  shown  in  Figure  2.  Separate  spectra  are 
shown  for  time  windows  about  the  Pg,  Lg,  and  Rg  phases,  respectively. 
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